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Numerical Simulations of a Protostellar Outflow Colliding with a Dense 

Molecular Cloud 

Chang Hyun Baek 1 ' 2 , Jongsoo Kim 3 , and Minho Choi 3 
ABSTRACT 

High-resolution SiO observations of the NGC 1333 IRAS 4A star-forming region 
showed a highly collimated outflow with a substantial deflection. The deflection was 
suggested to be caused by the interactions of the outflow and a dense cloud core. To 
investigate the deflection process of protostellar outflows, we have carried out three- 
dimensional hydrodynamic simulations of the collision of an outflow with a dense cloud. 
Assuming a power-law type density distribution of the obstructing cloud, the numer- 
ical experiments show that the deflection angle is mainly determined by the impact 
parameter and the density contrast between the outflow and the cloud. The deflection 
angle is, however, relatively insensitive to the velocity of the outflow. Using a numerical 
model with physical conditions that are particularly suitable for the IRAS 4A system, 
we produce a column-density image and a position-velocity diagram along the outflow, 
and they are consistent with the observations. Based on our numerical simulations, if 
we assume that the initial density and the velocity of the outflow are ~ 10 cm -3 and 
~ 100 km s _1 , the densities of the dense core and ambient medium in the IRAS 4A 
system are most likely to be ~ 10 5 cm~ 3 and ~ 10 2 cm~ 3 , respectively. We therefore 
demonstrate through numerical simulations that the directional variability of the IRAS 
4A outflow can be explained reasonably well using the collision model. 

Subject headings: ISM : jets and outflows — ISM : kinematics and dynamics — ISM : 
individual (NGC 1333 IRAS 4A) 

1. Introduction 

Radio interferometric observations in molecular lines show directional variability of jets in some 
star forming regions. The understanding of the variability of protostellar jets provides us with clues 
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to a variety of physical processes involved in the interaction of jets an d surrounding clouds. To 
explain the variabi li ty of protost ellar jet, several models were proposed (jEisloffel Sz Mundt 1 1 19971 ; 



Girart et al. Ill999i ; IChoi Il2001al ). These models can be classified into two categories. The first 



category includes external perturbation models, such as density gradients in the environment, flow 
instabilities (e.g., Kelvin-Helmholtz instability), side winds, and magnetic fields. The second cate- 
gory includes intrinsic variability models, such as the motion of the jet source perpendicular to the 
jet direction and the precession of the jet. Moreover, to understand the bending of the jet, espe- 
cially in the conte xt of young stellar objects (YSQs), numerical s imulations have been carried out 



by se veral authors (|de Gouveia Dal Pino Ill999l : iHurka et aZ.lll999l : iRaga Canto Ill995l ; iRaga et al 



2002). Three-dimensional smoothed particle hydr odynamic simulations o f the penetration of a jet 
into a dense stratified cloud were carried out by Ide Gouveia Dal Pino I (|1999i ) , and it was found 
that a radius ratio and a d ensity ratio between the jet and the cloud are important parameters in 
determining the outcome. IHurka et al\ (jl999l ) studied non-radiative jets bent by magnetic fields 
with a strong gradient. They showed that a low-velocity jet can be deflected by mag netic fields with 



a mod est strength, but a high-velocity jet needs very strong fields to be deflected. IRaga &: Canto 



()1995l ) presented analytical and two-dimensi onal numer i cal st udies of the collision of an Herbig- 
Halo jet with a dense molecular cloud core. IRaga et all (J2002J) studied the interaction between a 
radiative jet and a dense cloud through three-dimensional hydrodynamic simulations and generated 
Ha and H2I — 0S(1) emission maps. 

The NGC 1333 molecular cloud contains numerous YSOs and outflows and is a well-studied star 
formation region. A highly collimated outflow with a subst antial deflect ion angle was observed from 
high resolution SiO observations of NGC 1 333 IRAS 4A (IChoi l | 2005al). To invest i gate t he overall 
structure of the IRAS 4A outflow system, Ichoi I feoOll . l2005a! M and Ichoi etah feoOfil ) observed 
the outflow in t he C^O, 13 CO J = 1 -> 0, HCO+, HCN J = 1 -»• 0, SiO v = J = 1 -> 0, and H 2 
1 — S(l) lines. IChoi I (|2005al ) suggested that the sharp bend was caused by a collisi on between the 



northeastern outflow and a dense core in the ambient molecular cloud. And then, IChoi I (|2005b) 



confirmed that the obstructing cloud with a dense core is located just north of the bending point 
in images of the HCO + and HCN J = 1 — > lines. A sharp bending of the northeastern outflow 
in the NGC 1333 IRAS 4A system is unusual structure and provides a good example to study the 
bending mechanism of outflows in protostellar environment. 

In this study we have carried out three-dimensional hydrodynamic simulations on the inter- 
actions of an outflow with a dense molecular cloud surrounded by an ambient medium, in order 
to understand the deflection mechanism of the outflow and to construct a detailed model of de- 
flected outflows. We performed many simulations with different sets of physical parameters that 
can change the outflow direction, then picked up one model that gives a similar deflection angle 
with the observed one in the IRAS 4A system. Finally, we make a column-density image of shocked 
gas and a position-velocity diagram along the outflow. They are reasonably consistent with the 
observed data of the IRAS 4A northeastern outflow. In §2, we describe a numerical setup for 
the interaction of an outflow and a cloud. Simulation results are presented in §3, followed by the 
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summary and discussion in §4. 



2. Numerical Simulations 



We execute three-dimensional hydrodynamic simulations of outflow/cloud interacti ons with a 



grid- based hydrodynamic code based on the total variation diminishing (TVD) scheme (IRyu et al 
19931 ). We integrate numerically the hydrodynamic equations, 



dp „ 



(py) = o, 



dt 



(pv) + V • (pvv + pi) = 0, 



dE 

— + V.[(E + p)v] 



0, 



(1) 

(2) 
(3) 



where E = (1/2)/tv 2 +p/(7 — 1) and other notations have their usual meanings. For the adiabatic 
index, 7 = 5/3 is assumed. We did not take into account radiative cooling and heating processes, 
and chemical reactions, except for a test run explained below. The reason why we ignored them 
is that the deflection angl e, which is of our main intere st, is not sensitive to them. Two numerical 
simulations carried out by Ide Gouveia Dal Pino I (119991 ) showed that the deflection angle (6 ~ 40°) 
seen in the adiabatic simulation is slightly larger than the angle {0 ~ 35°) in the simulation 
with radiative co oling. To check their result , we also performed two simu lations with and without 



radia tive cooling ([Sutherland &: Dopitalll993l ; ISanchez-Salcedo et alj|2002l ) and heating (|Baek et al 



2005 ) processes. The structure of the outflow before and after collision in our two simulations are 
similar, except that a denser and shallower shell is formed at the head of the deflected outflow in 
the simulation with the radiative processes. The deflection angle shown in the adiabatic simulation 
is about y 3° larger than that shown in the simulation with the radiative processes. Both the 
results of Ide Gouveia Dal Pino I (|1999l ) and ours suggest that the radiative cooling does not play a 
major role in the outflow deflection. 

We setup a computational box whose minimum and maximum Cartesian coordinates are de- 
fined by -0.05 pc < x < 0, -0.025 pc < y < 0.025 pc, and -0.04 pc < z < 0.06 pc. The box is 
covered with 256 x 256 x 512 cells, so the numerical resolution in each direction becomes a constant, 
40 AU. A cloud with a radius R c = 0.02 pc is located at the origin of the coordinate system. An 
outflow with a sectional radius R Q = 300 AU is injected from the bottom x — y plane into the 
box with a direction parallel to the z-axis. Open boundary conditions are used on all boundaries, 
except for the bottom of x — y plane, on which a reflection condition is imposed. To explore the 
outflow/cloud interactions with different physical parameters, we vary the density distribution of 
a molecular cloud (an initially spherical molecular cloud with uniform density or power-law type 
density distribution), an impact parameter (from 3/10i? c to 10/10i? c ), and the density ratios among 
the outflow, cloud, and ambient medium. The initial ratios of the number density of outflow gas n a 
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to the number density of the ambient gas n a and the number density of outflow gas to the number 
density of cloud gas n c are defined as, respectively 

f] = — ,X=— ■ ( 4 ) 
n a n c 

Twelve simulations were done with different combinations of the model parameters, which are 
summarized in Table 1. 



3. Results 
3.1. Uniform Density Cloud Models 



To compare our simulation results with the results of iRaga et al\ ([2002), we performed nu- 
merical experiments on the interaction of an outflow with a uniform density cloud for different 
impact parameters, 4/10i? c , 6/10i? c , 8/10-R c , and 10/10i? c (See Table 1). The outflow with num- 
ber density n a = 10 cm -3 , temperature T a = 10 4 K, and velocity 100 km s _1 is injected toward an 
initially spherical cloud with uniform density n c = 10 3 cm -3 . The spherical cloud is surrounded by 
a homogeneous ambient medium of density n a = 1 cm -3 and temperature T a = 10 4 K. The cloud 
is in pressure equilibrium with the ambient medium. Figure 1 shows density and velocity fields on 
the y = plane for models UC01, UC02, UC03, and UC04. At t = 200 yr, the outflow has reached 
the cloud boundary and bores a hole through the cloud or gets deflected. In the case of the models 
with small impact parameters (UC01 and UC02) the incident outflow impinges on the surface of the 
cloud, makes a hole, and finally penetrates through the cloud. In the other case of the models with 
large impact parameters (UC03 and UC04), the incident outflow is deflected by the cloud, and then 
the deflected flow propagates with a wider opening angle. Depending on the impact parameter, the 
evolution of the interaction of an outflow with a uniform cloud can be classified into two cases as 
explained above. For the case of deflected outflow, the impact parameter is one of the important 
parameters that determines the deflection angle, which is the angle between the incident outflow 
direction and the deflected flow direction. Our results of the two uniform cloud models UC02 and 



UC03 are basically similar to those of model A of lRaga et al\ (|2002l ). 



3.2. Power-Law Type Cloud Models 

As we have pointed out in the introduction section, the main motivation of this paper is to 
make a numerical model of the deflected outflow seen in the NGC 1333 IRAS 4A system. So, 
whenever possible, we pick up physical quantities such as number density, and temperature from 
observations of the system. If not, we resort to the known quantities for the typical protostellar 
outflows and their environment. An outflow with number density n Q = 10 cnr 3 and temperature 
T Q = 10 4 K propagates through the uniform ambient medium with number density n a = 100 cm -3 
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and temperature T a = 10 K. The outflow interacts with an obstructing molecular cloud with a 
power-law type density profile, 

Ti 

nd= [l + A{r-r c W ^ r ^' ® 



where the distance from the cloud center r = \J x 2 + y 2 + z 2 , the cloud radius R c = 0.02 pc, 
the core radius r c = 400 AU, and the scaling factor A ~ 17 AU _1 . The number density and 
temperature in the center of the dense core are n c = 10 5 cm -3 and T c = 10 K, respectively. With 
the given number densities and temperatures for the outflow, cloud, and ambient medium, there 
is no pressure balance among those gas components. However, considering the fact that there is 
one order of magnitude difference between the largest sound speed with temperature 10 4 K and an 
outflow speed 100 km s" 1 , the motion induced by the pressure imbalance is negligible. 

Figure 2 depicts the density and velocity fields on the y = plane for models PC01, PC02, 
PC03, and PC04. For uniform density cloud models with smaller impact parameters (UC01 and 
UC02), the outflow makes a hole into and through the cloud (see Figure 1). For the power-law type 
density models, on the contrary, the outflow with a small impact parameter does not make a hole 
through the cloud but is deflected significantly. Furthermore, the deflection angle of the outflow 
increases as the impact parameter decreases. This result demonstrates that not only the impact 
parameter but also the density gradient of the cloud plays an important role in determining the 
deflection angle. We also measure deflection angles as a function of time from the simulation results 
of models PC01, PC02, PC03, and PC04, which is shown in Figure 3. The increase of deflection 
angles in the cases of small impact parameters (models PC01 and PC02) is quite significant, which 
is due to large density contrast between the outflow and the inner part of the obstructing cloud. It 
seems to us that, when the impact parameter is smaller than a critical value, the deflection angle 
increases with time and then keeps its angle after passing over the cloud. But outflow in models 
PC03 and PC04 is not efficiently disturbed by the obstructing cloud. For models PC03 and PC04 
the deflection angles are small, because the outflow interacts with the outer low density part of the 
cloud. After being deflected by the cloud, the outflow again keeps its small deflection angles for 
these models with large impact parameters. 

Figure 4 shows the time evolution of our fiducial model PC02, whose model parameters are 
n a = 10 cm -3 , rj = 0.1, and \ = 10 (see Table 1). At t = 700 yr, a bow shock generated by the 
outflow reaches the cloud boundary, and the velocity at the bow shock front is slightly changed due 
to the density gradient of the cloud boundary. The propagation direction of the outflow changes 
only slightly when the bow shock just passes through the cloud boundary. The larger deflection 
occurs at the nearest point from the center of the cloud at t ~ 1800 yr, and then the deflection 
angle tends to increase with time as the outflow passes over the cloud. At later times (t > 4200 yr), 
the deflected outflow moves into the ambient medium without any significant change of deflection 
angle. It is known that the o utflow velocities before and after the collision with the cloud has a 



relation (|Raga fc Canto Ill995l ). 

Vd ~ i>iCos6>, (6) 
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where is the velocity of the deflected outflow, Vi is the velocity of the incident outflow, and 6 is the 
deflection angle. At t = 4200 yr, a mean velocity of the deflected outflow vj, ~ 80 km s _1 measured 
from our simulation is approximately equal to the value Vd = 86.6 km s _1 (for = 100 km s -1 and 
6 ~ 30°) predicted by from the equation above. The deflection angle after t = 3800 yr for model 
PC02 is very similar to the observed value of the northeastern outflow of IRAS 4A. 

3.3. Models with Different Density Contrast and Outflow Velocity 

The effects of the outflow density and outflow velocity on the deflection angle can be examined 
by comparing the results from simulations with different outflow densities (models PC11 and PC12) 
and velocities (models PC21 and PC22). Models PC11 and PC12 have hundred (rj = 10, \ = 10~ 2 ) 
and ten (77 = 1, % = 10 -3 ) times larger outflow densities, respectively, than that of PC02 (77 = 0.1, 
X = 10" 4 ). And models PC21 and PC22 have two (200 km s" 1 ) and three (300 km s" 1 ) times 
larger outflow velocities, respectively, than that of PC02 (100 km s _1 ). Figure 5 shows the density 
and velocity fields for models PC11, PC12, PC21, and PC22. As shown in Figure 3, the light 
outflow, whose number density is smaller than the number density of ambient medium (77 < 0.1), 
experiences quite a large deflection, and its opening angle after collision is small. However the 
heavy outflow, whose number density is larger than or equal to the number density of ambient 
medium (77 > 1), experiences a slight deflection, and then it propagates with a wider opening angle 
(see Fig. 5 for models PC11 and PC12). The kinematic time scale, which is defined by the time 
required for the outflow to cross over the longer dimension of the computational box, is dependent 
upon the 77 parameter. In fact, the lower-left two panels of Figure 5 show that a bow shock of the 
model PC11 arrives the upper boundary at t = 2000 yr, whereas the bow shock of the model PC21 
does at t = 3500 yr. This is because, as the density of the outflow increases with respect to the 
density of the ambient medium, the outflow is less decelerated by the medium. The lower-right 
two panels show that the overall density and velocity structures of them are very similar to each 
other except for the different time epochs. The evolutionary time of the faster outflow model PC22 
should naturally be shorter than that of model PC21. This result show that the outflow velocity 
does not influence the deflection angle. Considering the three observational facts of the IRAS 4A 
outflow, (i) the deflection angle of ~ 30°, (ii) the well-collimated flow, and (iii) the kinematic time 
scale of the outflow, PC02 is the fiducial model that describes the IRAS 4A system most closely. 

3.4. Comparison between Observational and Simulated Results 

The deflection angle of the outflow in model PC02 may be compared directly with the observed 
angle of the IRAS 4A SiO outflow (Figure 6). To make the simulated outflow image shown in the 
left box of Figure 6, we first pick up shocked gas whose temperature is higher than 3 times the 
initial outflow temperature (T > 3T ), then integrate the shocked gas along the y-axis. We can 
find, in the simulated outflow image, the well collimated flow and the enhanced column density 
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at the bending point and the terminal part of the outflow, which are very similar to the observed 
SiO image. Moreover, the time epoch of the simulated outflow image is similar to the kinem atic 
timescale of IRAS 4A outflow determined by the SiO main outflow observation (jChoi Il2005al ). 



The origin of the SiO molecules is not well-known. They can be either injected into the flow 
at the base of the outflow near the protostar or gradually mixed in from the ambient cloud. To 
trace the outflow and cloud separately, we have implemented a method of the Largrangian tracer 
(jJones et aZ.1 119961 ) into our Eulerian hydrodynamic code. From this Largrangian tracer, we find 
that the shocked gas mainly consists of material from the outflow, although a small amount of 
material from the cloud and ambient medium is entrained in the deflected shock gas. From the 
information on the amount of shocked gas, we infer that most of the observed SiO gas might be 
from the outflow. However, our statement needs to be confirmed by simulations that explicitly 
includes the chemistry of SiO and related species. 

The complicated kinematics of the deflected flow was revealed using position-velocity diagram. 
Figure 7 depicts the column density image and the position-velocity diagram along the solid line 
in the column density image calculated from model PC02. When we make the position- velocity 
diagram, we consider that the y = plane coincides with the plane of the sky and take into account 
the inclination angle 9° of the outflow with respective to the plane of the sky, which was chosen from 
two assumptions ; 1) the typical velocity of the outflow is 100 km s _1 , 2) the velocity interval in the 
position-velocity diagram is ~ 20 km s^ 1 . The position-velocity diagram shows that the line-of- 
slight velocity of the deflected part is slower and more turbulent than that of the undeflected part. 
This implies that the outflow is slowed down and heated by the collision with a dense molecular 
core and the interaction with the ambient medium. The broad-velocity structure near the head of 
the deflected o utflow is similar to that of the position-velocit y diagram of the SiO observation. It 
was suggested (|Mac Low h Klessen 1 12004 ) and investigated (jNakamura Li 1120071 ) that outflows 
might be one of the mechanisms that drive turbulence in molecular clouds. In our simulations, even 
though they are single collision experiments of the outflow and the cloud, we demonstrate that a 
turbulent flow with a velocity dispersion of the order of 10 km sec -1 can be easily generated. 



4. Summary and Discussion 

We have performed three-dimensional hydrodynamics simulations to study the interaction of 
an outflow with a dense molecular cloud. Since the most evident information from observations is 
the deflection angle, we focus on the deflection angle in the simulations by varying four parameters: 
the impact parameter, the outflow density contrast with the ambient medium (77), the contrast 
between outflow and cloud density (x), and the velocity of outflow. The main results can be 
summarized as follows. 

• When the impact parameter is smaller than 0.6 times the cloud radius, the outflow bores a 
hole into the uniform cloud. On the contrary, the outflow with the impact parameter as small 
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as 0.3 is deflected by the cloud with the power-law type density distribution. 

When the obstructing cloud has a power-law density distribution, the deflection angle of the 
outflow is determined by the impact parameter and the density contrast between the outflow 
and the cloud (x). 

The deflection angle is relatively insensitive to the velocity of outflow. 

The initially light outflow (77 < 1) can be more collimated than the heavy one (77 > 1). 



Choi I (|2005al ) suggested that the NGC 1333 IRAS 4A outflow was deflected as a result of a 
collision with a dense cloud core and provided several lines of evidence. They are (1) the asymmetric 
morphology of the bipolar outflow, (2) the good collimation and complicated kinematics of the 
deflected flow, (3) the low-velocity emission from the molecular gas near the bend, and (4) the 
enhancement of SiO emission in the deflected flow. In this study, we confirmed these evidences 
through numerical experiments. Especially, the deflection angle (~ 30°) of collimated outflow and 
the enhancement and complicated kinematics of shocked gas in the deflected region produced from 
our numerical simulations are very similar to those seen in the SiO observations of the IRAS 4A 
region. If we take model PC02 as the most appropriate numerical model for the IRAS 4A outflow, 
we can infer that this outflow was significantly deflected by the dense molecular core ~ 2000 years 
ago and the ratio of the outflow density to the density at the cloud center is very low ~ 10 . 
Moreover, if we assume that the initial density and the velocity of the outflow are ~ 10 cm~ 3 and 
~ 100 km s , the densities of the dense core and ambient medium in the IRAS 4A system are 
most likely to be ~ 10 5 cm -3 and ~ 10 2 cm -3 , respectively. 

Although our simulations do not consider the radiative cooling and chemical reactions, our 
results provide some insights to the interactions of an outflow with a dense cloud. To make a more 
sophisticated numerical model compatible with the SiO observations, we need to include chemistry, 
especially related to the SiO molecule. In the near future, we will execute numerical simulations 
which include the radiative cooling and chemical reactions related to the SiO molecule and then 
try to directly compare the simulated SiO emission image with the observed one. 
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Table 1. Model Parameters for Simulations 



Model 


Impact Parameter 


Cloud Density Distribution 


„ _ no 
1 n a 


A/ n ° ^ 


v b 


UC01 


4/10 R c 


uniform 


10 


io- 2 


100 


UC02 


6/10 R c 


uniform 


10 


io- 2 


100 


UC03 


8/10 R c 


uniform 


10 


IO" 2 


100 


UC04 


10/10 R c 


uniform 


10 


io- 2 


100 


PC01 


3/10 R c 


power-law type 


0.1 


io- 4 


100 


PC02 


4/10 R c 


power-law type 


0.1 


io- 4 


100 


PC03 


5/10 R c 


power-law type 


0.1 


10~ 4 


100 


PC04 


6/10 R c 


power-law type 


0.1 


io~ 4 


100 


PC11 


4/10 R c 


power-law type 


10 


io- 2 


100 


PC12 


4/10 R c 


power-law type 


1 


IO" 3 


100 


PC21 


4/10 R c 


power-law type 


0.1 


io- 4 


200 


PC22 


4/10 R c 


power-law type 


0.1 


io- 4 


300 



a For the power-law type density models, n c represents the central number density of 
a cloud. 

b The outflow velocity in units of km s _1 . 
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Fig. 1. — Density and velocity fields of uniform density cloud models UC01, UC02, UC03, and 
UC04. Gray images for density fields and arrows for velocity fields on the y = plane (which 
includes the outflow axis and the center of the spherical cloud) are shown in the upper panels at 
a time epoch t = 1000 yr and the lower ones at t = 3000 yr. The gray scale bar is in logarithmic 
units. The longest arrow corresponds to 100 km s . 
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Fig. 2. — Density and velocity fields of power-law type density cloud models PC01, PC02, PC03, 
and PC04. The model name is given at top of each upper panel. The upper and lower panels 
show gray images for density fields and arrows for velocity fields at time epochs 1000 yr and 4500 
yr, respectively. The gray scale bar is in logarithmic units. The longest arrow corresponds to 
100 km s _1 . 
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Fig. 3. — Evolution of deflection angle for models PC01 (dotted line; filled circle), PC02 (solid line; 
open circle), PC03 (dashed line; filled box), and PC04 (long-dashed line; filled triangle). 
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Fig. 4. — The time evolution of density and velocity fields for model PC02. Gray images for density 
fields and arrows for velocity fields on the y = plane are shown. A two digit number in each panel 
shows a time epoch in units of 100 yr. 



-15- 




Fig. 5. — Density and velocity fields for models PC11, PC12, PC21, and PC22. Gray images for 
density fields and arrows for velocity fields on the y = plane are shown. A gray scale bar in 
logarithmic units is included in the far left side. The longest arrow corresponds to a speed around 
300 km s _1 . A two digit number in each panel shows a time epoch in units of 100 yr. 



Fig. 6. — Comparison of observed and simulated outflows (from lBaek et aZ.ll2007l ) . The right box 
contains a SiO image toward the northeastern region of NGC 1333 IRAS 4A. The left box shows 
an image of the column density of shocked gas at a time epoch of 4200 yr of the numerical model 
PC02. The deflection angles seen in both boxes are very similar. 
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Fig. 7. — The column-density image and position- velocity diagram of model PC02 in Table 1. The 
position-velocity diagram in the right panel is plotted along the solid line in the left panel. The 
intersection of the solid and dashed lines in the left panel corresponds to the bending point of the 
outflow. 



